Theory of orbital magnetoelectric response 

Andrei Malashevich 1 , Ivo Souza 1 , Sinisa Coh 2 and David 
Vanderbilt 2 

1 Department of Physics, University of California, Berkeley, CA 94720-7300, USA 

2 Department of Physics & Astronomy, Rutgers University, Piscataway, NJ 
08854-8019, USA 

E-mail: andreim@berkeley.edu 

Abstract. We extend the recently-developed theory of bulk orbital magnetization 
to finite electric fields, and use it to calculate the orbital magnetoelectric response 
of periodic insulators. Working in the independent-particle framework, we find that 
the finite-field orbital magnetization can be written as a sum of three gauge-invariant 
contributions, one of which has no counterpart at zero field. The extra contribution 
is collinear with and explicitly dependent on the electric field. The expression for 
the orbital magnetization is suitable for first-principles implementations, allowing 
to calculate the magnetoelectric response coefficients by numerical differentiation. 
Alternatively, perturbation-theory techniques may be used, and for that purpose we 
derive an expression directly for the linear magnetoelectric tensor by taking the first 
field-derivative analytically. Two types of terms are obtained. One, the 'Chern-Simons' 
term, depends only on the unperturbed occupied orbitals and is purely isotropic. The 
other, 'Kubo' terms, involve the first-order change in the orbitals and give isotropic 
as well as anisotropic contributions to the response. In ordinary magnetoelectric 
insulators all terms are generally present, while in strong Zi topological insulators only 
the Chern-Simons term is allowed, and is quantized. In order to validate the theory 
we have calculated under periodic boundary conditions the linear magnetoelectric 
susceptibility for a 3-D tight-binding model of an ordinary magnetoelectric insulator, 
using both the finite-field and perturbation-theory expressions. The results are in 
excellent agreement with calculations on bounded samples. 
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1. Introduction 

In insulating materials in which both spatial inversion and time-reversal symmetries 
are broken, a magnetic field B can induce a first-order electric polarization P, and 
conversely an electric field £ can induce a first-order magnetization M [1, 2]. This 
linear magnetoelectric (ME) effect is described by the susceptibility tensor 



dP d 

Otda = 



0B a 



dM a 



(1) 

£=0 



B=0 

where indices label spatial directions. This tensor can be divided into a "frozen-ion" 
contribution that occurs even when the ionic coordinates are fixed, and a "lattice- 
mediated" contribution corresponding to the remainder. Each of these two contributions 
can be decomposed further according to whether the magnetic interaction is associated 
with spins or orbital currents, giving four contributions to a in total. 

All of those contributions, except the frozen-ion orbital one, are relatively 
straightforward to evaluate, at least in principle, and ab initio calculations have 
started to appear. For example, the lattice-mediated spin-magnetization response was 
calculated in [3] for Cr 2 3 and in [4] for BiFe0 3 (including the strain deformation 
effects that are present in the latter), and calculations based on the converse approach 
(polarization response to a Zeeman field) were recently reported [5]. One generally 
expects the lattice-mediated couplings to be larger than the frozen-ion ones, and insofar 
as the spin-orbit interaction can be treated perturbatively, interactions involving spin 
magnetization are typically larger than the orbital ones. However, we shall see that there 
are situations in which the spin-orbit interaction cannot be treated perturbatively, and 
in which the frozen- ion orbital contribution is expected to be dominant. Therefore, it 
is desirable to have a complete description which accounts for all four contributions. 

The frozen-ion orbital contribution is, in fact, the one part of the ME susceptibility 
for which there is at present no satisfactory theoretical or computational framework, 
although some progress towards that goal was made in two recent works [6, 7]. Following 
Essin et al. [7] we refer to it as the "orbital magnetoelectric polarizability" (OMP). For 
the remainder of this paper, we will focus exclusively on this contribution to (1), and 
shall denote it simply by a. Accordingly, the symbol M will be used henceforth for the 
orbital component of the magnetization. 

The question we pose to ourselves is the following: what is the quantum-mechanical 
expression for the tensor a of a generic three-dimensional band insulator? We note that 
the conventional perturbation-theory expression for a [8, 9] does not apply to Bloch 
electrons, as it involves matrix elements of unbounded operators. The proper expressions 
for P [10] and M [11, 12, 13, 14] in periodic crystals have been derived, but so far only 
at B = and £ = respectively. The evaluation of equation (1) remains therefore an 
open problem. 

Phenomenologically, the most general form of a is a 3 x 3 matrix where all nine 
components are independent. Dividing it into traceless and isotropic parts, the latter is 
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conveniently expressed in terms of a single dimensionless parameter 9 as 

The presence of an isotropic ME coupling is equivalent to the addition of a term 
proportional to 9£ ■ B to the electromagnetic Lagrangian. Such a term describes "axion 
electrodynamics" [15] and (2) may therefore also be referred to as the "axion OMP." 
The electrodynamic effects of the axion field are elusive (in fact, the very existence of 
a 6 was debated until recently: see [16, 17] and references therein). For example, in a 
finite, static sample cut from a uniform ME medium those effects are only felt at the 
surface [15, 18]. In particular, a e gives rise to a surface Hall effect [19]. 

An essential feature of the axion theory is that a change of 9 by In leaves the 
electrodynamics invariant [15]. The profound implications for the ME response of 
materials were recognized by Qi et al. [6] , and discussed further by Essin et al. [7] . These 
authors showed that there is a part of the isotropic OMP which remains ambiguous up 
to integer multiples of 2ir in the corresponding 9 until the surface termination of the 
sample is specified. For example, a change by 2nn occurs if the surface is modified by 
adsorbing a quantum anomalous Hall layer. Hence this particular contribution to 9 can 
be formulated as a bulk quantity only modulo a quantum of indeterminacy, in much 
the same way as the electric polarization P [10, 20]. A microscopic expression for it 
was derived in the framework of single-particle band theory by the above authors. It 
is given by the Brillouin-zone integral of the Chern-Simons form [21] in /c-space, which 
is a multivalued global geometric invariant reminiscent of the Berry-phase expression 
for P [10]. We denote henceforth this "geometric" contribution to the OMP as the 
Chern-Simons OMP (CSOMP). 

A remarkable outcome of this analysis is the prediction [6] of a purely isotropic 
"topological ME effect," associated with the CSOMP, in a newly-discovered class of 
time-reversal invariant insulators known as Z 2 topological insulators [22, 23, 24]. As a 
result of the multivaluedness of 9, the presence of time-reversal symmetry in the bulk, 
which takes 9 into —9, is consistent with two solutions: 9 = 0, corresponding to ordinary 
insulators, and 9 = it, corresponding to strong Z 2 topological insulators.^ The latter 
case is non-perturbative in the spin-orbit interaction, and 9 = n amounts to a rather 
large ME susceptibility (in Gaussian units it is 1/4-7T times the fine structure constant, 
or ~6xl0~ 4 , to be compared with ~lxl0 -4 for the total ME response of &2O3 at low 
temperature [25]). 

It is not clear from these recent works, however, whether the isotropic CSOMP 
constitutes the full OMP response of a generic insulator. It does appear to do so for the 
tight-binding model studied in [7] , whose ME response was correctly reproduced by the 
Chern-Simons expression even when the parameters were tuned to break time-reversal 

X An analogous situation occurs in the theory of polarization: inversion symmetry, which takes P into 
P, allows for a nontrivial solution which does not include P = in the "lattice" of values [20]. An 
important difference is that while 9 is a directly measurable response, only changes in P are detectable, 
so that the experimental implications of the nontrivial solution are less clear in this case. 
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and inversion symmetries (i.e., for generic not equal to or ir). On the other hand, 
other considerations seem to demand additional contributions. For example, it is not 
difficult to construct tight-binding models of molecular crystals in which it is clear that 
the OMP cannot be purely isotropic. 

In this work we derive, using rigorous quantum-mechanical arguments, an 
expression for the OMP tensor a of band insulators, written solely in terms of bulk 
quantities (the periodic Hamiltonian and ground state Bloch wavefunctions, and their 
first-order change in an electric field). We restrict our derivation to non- interacting 
Hamiltonians, as the essential physics we wish to describe occurs already at the 
single-particle level. We find that in crystals with broken time-reversal and inversion 
symmetries there are, in addition to the CSOMP term discussed in [6, 7], extra terms 
which generally contribute to both the trace and the traceless parts of a. 

Our theoretical approach closely mimics one type of ME response experiment: a 
finite electric field £ is applied to a bounded sample, and the (orbital) magnetization 
is calculated in the presence of the field. Then the thermodynamic limit is taken at 
fixed field. This key step in the derivation must be done carefully, so that crucial 
surface contributions are not lost in the process, and here we follow the Wannier-based 
approach of references [12, 13], adapted to £ ^ 0. Finally the linear response coefficient 
otda — dM a /dSd is extracted in the limit that £ goes to zero. 

In a concurrent work by Essin, Turner, Moore, and one of us [26] an alternative 
approach was taken, which is closer in spirit to the calculation in [10] of the change 
in polarization as an integrated current: the adiabatic current induced in an infinite 
crystal by a change in its Hamiltonian in the presence of a magnetic field is computed, 
and then expressed as a total time derivative. The two approaches are complementary 
and lead to the same expression for a, illuminating it from different angles. 

The paper is organized as follows. In section 2 we derive the bulk expression for 
M(£), and reorganize it into three gauge- invariant contributions, one of which yields 
directly the CSOMP response. The gauge-invariant decomposition of M(£) is done 
at first in fc-space for periodic crystals, and then also for bounded samples working in 
real space. In section 3 we derive a fc-space formula for the OMP tensor a by taking 
analytically the field-derivative of M{£). Numerical tests on a tight-binding model of 
a ME insulator are presented at appropriate places throughout the paper in order to 
validate the bulk expressions for M{£) and a. In Appendix A we describe the tight- 
binding model, as well as technical details on how the various formulas are implemented 
on a Appoint grid. Appendix B and Appendix C contain derivations of certain results 
given in the main text. 
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2. Orbital magnetization in finite electric field 

2.1. Preliminaries 

The orbital magnetization M is defined as the orbital moment per unit volume, 



Here e > is the magnitude of the electron charge, V is the sample volume, and 
are the occupied eigenstates. While this expression can be directly implemented 
when using open boundary conditions, the electronic structure of crystals is more 
conveniently calculated and interpreted using periodic boundary conditions, in order 
to take advantage of Bloch's theorem. This poses however serious difficulties in dealing 
with the circulation operator r x v, because of the unbounded and nonperiodic nature 
of the position operator r. These subtle issues were fully resolved only recently, with 
the derivation of a bulk expression for M directly in terms of the extended Bloch 
states [11, 12, 13, 14]. 

In previous derivations the crystal was taken to be under shorted electrical boundary 
conditions. We shall extend the derivation given in [12, 13] to the case where a static 
homogeneous electric field £ is present, so that the full Hamiltonian reads 



The derivation, carried out for an insulator with N valence bands within the 
independent-particle approximation, involves transforming the set of occupied 
eigenstates of % into a set of Wannier-type (i.e., localized and orthonormal) orbitals 
\wj) and expressing M{£) in the Wannier representation. This is done at first for a 
finite sample cut from a periodic crystal, and eventually the thermodynamic limit is 
taken at fixed field. 

Before continuing, two remarks are in order. First, the assumption that it is possible 
to construct well-localized Wannier functions (WFs) spanning the valence bands is only 
valid if the Chern invariants of the valence bandstructure vanish identically [27]. This 
requirement is satisfied by normal band insulators as well as by Z2 topological insulators, 
but not by quantum anomalous Hall insulators [28], which thus far remain hypothetical. 
Second, because of Zener tunnelling, an insulating crystal does not have a well-defined 
ground state in a finite electric field. Nonetheless, upon slowly ramping up the field to 
the desired value, the electron system remains in a quasistationary state which is, for 
all practical purposes, indistinguishable from a truly stationary state. This is the state 
we shall consider in the ensuing derivation. As discussed in [29, 30], it is Wannier- and 
Bloch-representable, even though the Hamiltonian (4) is not lattice-periodic. 

2.2. k-space expression 

Our derivation of a A;-space (bulk) expression for M{£) is carried out mostly in real 
space, using a Wannier representation. It is only in the last step that we switch to 




(3) 



H = H u + e£-r. 



(4) 
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reciprocal space, by expressing the crystalline WFs \Rn) in terms of the cell-periodic 
Bloch functions \u n k) via [31] 

\Rn) = V c J[dk]e ik < r -V\u nk ), (5) 

where R is a lattice vector, V c is the unit-cell volume, [dk] = d 3 /c/(27r) 3 , and the integral 
is over the first Brillouin zone. 

We begin with a finite sample immersed in a field £, divide it up into an interior 
region and a surface region, and assign each WF to either one. The boundary between 
the two regions is chosen in such a way that the fractional volume of the surface 
region goes to zero as V — > oo, but deep enough that WFs near the boundary are 
bulk-like. Following [12, 13], equation (3) for the orbital magnetization can then be 
rewritten as an interior contribution plus a surface contribution, denoted respectively 
as the "local circulation" (LC) and the "itinerant circulation" (IC). Remarkably, in 
the thermodynamic limit both can be expressed solely in terms of the interior-region 
crystalline WFs, or equivalently, in terms of the bulk Bloch functions, as shown in the 
above references at £ = and below for £ ^ 0. Specifically, we shall show that 

M = M LC + M IC '° + M 1C ' £ , (6 a) 

where 

N 

M a LC = - 7 e a6c Im^ [dk](d b u nk \H° k \d c u nk ) (66) 

n 

is the contribution from the interior WFs, 

N 

Ml C >° = -ie abc lm I [dfc] (d b u nk \d c u mk ) H° mnk (6c) 

nm 

is the part of the surface contribution coming from the zero-field Hamiltonian, and 

N 

Ml°> £ = -7e abc Im^ / [dk}(d b u nk \d c u mk )e£ ■ A mnk (6<f) 

nm 

is the part of the surface contribution coming from the electric field term in the 
Hamiltonian (4). In the above expressions 7 = —e/(2hc), 

Hi = e~ ikr U°e ikr , (7) 

H mnk = (umk\H k \u nk ) , (8) 

and A mnk is the Berry connection matrix defined in equation (14) below. 

Having stated the result we now present the derivation, starting with the interior 
contribution M . Using [rj,rj] = 0, the velocity operator v = (i/h)[H, r] becomes 
(i/K)[H°, r], so that the circulation operator r x v is unaffected by the electric field. It 
immediately follows that the local circulation part M LC is given in terms of the field- 
polarized states \u nk ) by the same expression, equation (66), as was derived in [13] for 
the zero-field case. 
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Consider now the contribution M lc = M 1C '° + M 1C ' £ from the surface WFs \w s ). 
For large samples it takes the form [13] 

surf 

M IC = --4 7T7 y r s x v s , (9) 
2cN c V c ^ y ' 

where N c is the number of crystal cells of volume V c , r s = (w s \r\w s ), and 

2 

v s = (w s \v\w s ) = -Im(w s \r'H\w s ). (10) 

Note that %\w s ) already belongs to the occupied manifold spanned by P = 
^° cc \wj)(wj\, since we assume a (quasi)stationary state. Thus we can insert a P 
between r and T-L above, and using (4) we obtain 

N 

^-Elii+i))- ( n ) 

3 

where v ^ = (2/h)lm[r S j / H^ s \ is the same as in [12, 13] and v^^ = (2e/h)lm[r S j(rj S - £)] 
is a new term. 

The reasoning [12, 13] by which M IC can be recast in terms of the bulk WFs 
\Rn) relies on the exponential localization of the WFs and on certain properties of 
(antisymmetry under j -h- s and invariance under lattice translations deep inside the 
crystallite) which are shared by v^ s y Hence we can follow similar steps as in those 
works, arriving at 

N 

M l°' £ = E E vfo m ,Rn } , b Rc (12) 

c R mn 

and similarly for M^ c '° with v° substituting for v £ . The latter is identical to the 
expression for M^ c valid at £ = [12, 13], and upon converting to /c-space becomes 
(6c). ' 

Let us now turn to M^ c ' £ and write (12) as (e 2 /2chV c )e abc £ d lmW bdtC where 

N 

W bd ,c = ^^(Rn\r b \0m)(0m\r d \Rn)R c . (13) 

R ran 

In order to recast this expression as a A;-space integral it is useful to introduce the N xN 
Berry connection matrix 

A m nk,b = i(u m k\dbU n k) = A^ mkb , (14) 

where d b = d/dk b . It satisfies the relation [31, 32] 

(Rn\r b \ Om) = V C J [dk]A nmk:b e ik - R . (15) 



We also need 



R c (Rn\r d \0m) = iV c J [dk](d c A nmk4 )e ik - R , (16) 
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Figure 1. The zz and zy components of the OMP tensor a of the tight-binding 
model described in Appendix A, as a function of the parameter ip. The two lower 
bands are treated as occupied. Solid line: extrapolation from finite-size samples using 
numerical differentiation of the finite- field magnetization calculated from (3). Dashed 
line: numerical differentiation of the finite- field magnetization calculated using (66)- 
(6rf) discretized on a /c-space grid. Open circles: linear-response calculation in fc-space 
using discretized versions of (47a)-(47c). 



and we arrive at (6(f). 

The sum of equations (66)— (6cQ gives the desired A;-space expression for M(£). In 
the limit S — > the term M 1C ' £ vanishes, and equation (31) of [13] is recovered. 

We have implemented (66)— (6<i) for the tight-binding model of Appendix A. Since 
for small electric fields M{£) differs only slightly from M(0), in order to observe 
the effect of the electric field we consider differences in magnetization rather than the 
absolute magnetization. Therefore, in all our numerical tests we evaluated the OMP 
tensor a da - With the help of (66)-(6cf) we calculated it as AM„/A^, using small fields 
£ d = ±0.01. We then repeated the calculation on finite samples cut from the bulk 
crystal, using (3) in place of (66)— (6cQ. Figure 1 shows the value of the zz and zy 
components of a plotted as a function of the parameter <p, the phase of one of the 
complex hopping amplitudes (see Appendix A for details). The very precise agreement 
between the solid and dashed lines confirms the correctness of the /c-space formula. The 
same level of agreement was found for the other components of a. 



which follows from (15). Using these two relations, (13) becomes 




mn 



(17) 



2.3. Gauge-invariant decomposition 
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2.3.1. Periodic crystals Equation (6a) for M{£) is valid in an arbitrary gauge, that 
is, the sum of its three terms given by (6o)-(6d) - but not each term individually - 
remains invariant under a unitary transformation 

N 

\u„ k ) -> ^2 \ U mk)U m nk (18) 
m 

among the valence-band states at each k. In order to make the gauge invariance of (6a) 
manifest, it is convenient to first manipulate it into a different form, given in terms of 
certain canonical objects which we now define. We begin by introducing the covariant 
/c-derivative of a valence state [30], 

\d b u n k) = Qk \d b u nk ), (19) 

where Qk = 1 — Pk and 

N 

Pk = ^2\ujk)(u jk \- (20) 

3=1 

The covariant and ordinary derivatives are related by 

N 

\d b u nk ) = \d b u nk ) - A mnkjb \u m k) . (21) 

m 

The generalized metric-curvature tensor is [31] 

Fnmk,hc = (p b U n k\d c U m k) — F m nk,c b ' (22) 

Viewed as an N x iV matrix over the band indices, F is gauge-covariant, changing as 
Fnmk, b c -+ (UlF ktbc U k ) (23) 

under the transformation (18). We also note the relation 

{d b u n k\d c u m k) F nm k )b c ~\~ (^4fc,fe^4fc,c)n?Tf (24) 

We shall make use of two more gauge-covariant objects, 

H nmk,b = i{Unk\Hl\d b U mk ) (25) 

and 

H nmk, bc = (d b u nk \H^\d c u mk ) , (26) 
which enter the relation 

(d b u nk \H° k \d c u mk ) = H° nmkibc + \A k , b Hl c + (Hl b y A k , c + A k , b H° k A k , c ] . (27) 

L J nm 

Coming back to equations (6a)-(6d), for M^ c we use (27) and for M^ c we use (24), 
leading to 

M a = - 7 e a6c J [dk] Im tr [#° + 2A b H° c + H°F bc + e£ d A d F bc + eS d A d A b A c j , (28) 
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where "tr" denotes the electronic trace over the occupied valence bands and we have 
dropped the subscript k. The second term can be rewritten using 

Hnm,c = ~ e ^dFnm,dc- (29) 

(To obtain this relation start from the generalized Schrodinger equation satisfied by the 
valence states at £ ^ [33], 

N 

H°\u n ) = ( H mn + e£ • A mn )\u m ) - ie£ d \d d u n ), (30) 

m 

and multiply through by (d c u m \.) 
Let us define the quantities 



M a 



LC 



M: 



ic 



-je abc J[dk}lmti[H bc ], (31) 
-ie abc J[dk]lmti[H°F bc }, (32) 



and 



M a cs = -e-fe abc £ d J [dk]lm tr[2 A b F cd + F bc A d + A b A c A d }. (33) 

The total magnetization is given by their sum 

M a = M^c + mic + M cs_ ( 34 ^ 

Referring to (22) and (26) the first two terms read, in a more conventional notation, 
_ r N 

M^ C = -ie abc / [dk] Im(d b u nk \H° k \d c u nk ) (346) 

J n 

and 

r N 

M a ° = - 7 e abc / [dk] Im (^mfc^mfcl^nfc)) ■ (34 c) 

nm 

These are the only terms that remain in the limit £ — > 0, in agreement with equation (43) 
of [13]. At finite field they depend on £ implicitly via the wavefunctions. 

We now show that the term M cs , which gathers all the contributions with an 
explicit dependence on £, can be recast as 



M a cs = e 1 £ a J [dk]t ijk ti 



AidjA^ —AiAjAj- 



(34 d) 



To do so it is convenient to introduce the Berry curvature tensor 

Qnm,ab iFnm,ab ^Fnm,ba ^nrn^bai (35) 

where F nm>ab was defined in (22). A few lines of algebra show that 

^nm,ab d a A nmb d b A nma i[_A a , A b ] nm . (36) 

In order to go from (33) to (34(f), use (35) to write Im tr[i^ c AJ as — | trLA^bc] and 
— 21m tr[A fe F rfc ] as trLA b fi dc ], and then replace VL nm>bc in these expressions with e abc Q nm:a , 
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where fl nm , a — \ e abc^nm,bc 1S the Berry curvature tensor written in axial- vector form. 
This leads to 

M a cs = e 7 J [dfc](£ a tr[fi • A] - S d e abc lm tr[A b A c A d ]). (37) 

The first term is parallel to the field, and can be rewritten with the help of (36): 

tr[fi • A] = e ijk tT[AidjA k - \AiAjA k \. (38) 

While not immediately apparent, the second term in (37) also points along the field. To 
see this, write 

^2 £de a6c Im tr[A b A c A d ] = £ a ^ e abc lm tr[A a A b A c ] + e « fecIm tr[A b A c A d ], (39) 

bed be d^a be 

where we suspended momentarily the implied summation convention. The last term 
vanishes because the factor e abc forces d ^ a to equal either b or c, producing 
terms such as Im tr^A^Ac] which vanish identically as A b is Hermitian. Rewriting 
£ a Y.bc e a6clm trL4 a A b A c ] as (£ a /3) J2ijk e ijk lm trlAiAjAk] and restoring the summation 
convention, we arrive at (34 cf) . 

Equations (346)-(34d), which constitute the main result of this section, are 
separately gauge- invariant. For M LC and M 1C this is apparent already from (31) and 
(32), whose integrands are gauge-invariant, being traces over gauge-covariant matrices. 
In contrast, equation (34c?) for M cs only becomes invariant after taking the integral 
on the right-hand-side over the entire Brillouin zone (the integrand being familiar from 
differential geometry as the Chern-Simons 3-form [34, 21]). 

The Chern-Simons contribution (34(f) has several remarkable features: (i) as already 
noted, it is perfectly isotropic, remaining parallel to £ for arbitrary orientations of £ 
relative to the crystal axes; (ii) being isotropic, it vanishes in less than three dimensions, 
which intuitively can be understood because already in two dimensions polarization 
must be in the plane of the system and magnetization must be out of the plane; (iii) for 
N > 1 valence bands it is a multivalued bulk quantity with a quantum of arbitrariness 
(e 2 /hc)£ a , a fact that is connected with the possibility of a cyclic adiabatic evolution 
that would change (47a) below for 6 by 2n [6]. 

We have repeated the calculation of the OMP presented in figure 1 using (346)- 
(34(f) instead of (66)-(6d), finding excellent agreement between them. The electric field 
derivative of the decomposition (34 a) gives the corresponding decomposition of the OMP 
tensor (1), 

a = 5 LC + 5 IC + a cs , (40) 

where each term is also gauge-invariant. The zz components of these terms are plotted 
separately in figure 2. 

2.3.2. Finite samples It is natural to ask whether the gauge-invariant decomposition of 
the orbital magnetization given in equation (34a) can be made already for finite samples, 
before taking the thermodynamic limit and switching to periodic boundary conditions. 
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Figure 2. Decomposition of the a zz curve in figure 1 into the gauge- invariant 
contributions (solid lines), a l p z (dashed line), and a^f (dotted line), calculated in 
/c-space using finite differences in 5. Symbols denote the same contributions evaluated 
for bounded samples, also using finite differences. 



This has previously been done in the case £ = 0, where M cs = and M LC and M 1C 
take the form [35] 

M\ c = ^^ abc lm Tr [Pr b Qn°Qr c } (41a) 

and 

Ml c = ^e a6c ImTr [PU°Pr b Qr L ]. (416) 

Here P and Q — 1 — P are the projection operators onto the occupied and empty 
subspaces, respectively, and "Tr" denotes the electronic trace over the entire Hilbert 
space. These two expressions, which are manifestly gauge-invariant, remain valid at 
finite field, reducing to (346) and (34c) in the thermodynamic limit. 

We now complete this picture for 8 ^ by showing that the remaining contribution 
M cs = M — M LC - M lc can also be written in trace form, as 

M a cs = -^S a e ijk lm Tr [Pr.Pr.Pn]. (41c) 
We first recast the orbital magnetization (3) as 

M a = -^ytabcTr [Pr b v c ] = ^L_e abc lm Tr [Pr b V°r c ] (42) 
and then subtract (41a) and (416) from it to find, after some manipulations, 

Replacing VP with H — eE^d and using QHP = 0, 

e 2 

M a cs = -— e^Im Tr [Pr d Pr b Pr c ] . (44) 



M a cs = -^e a6c ImTr [Qn°Pr b Pr c ]. (43) 
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The imaginary part of the trace vanishes if any two of the indices 6, c, or d are the 
same, and therefore d must be equal to a. Using the cyclic property we conclude that 
all non- vanishing terms in the sum over b and c are identical, leading to (41c). This 
part of the field-induced magnetization is clearly isotropic, with a coupling strength (see 
equation (2)) given by 

Atx 2 

e cs = -^7^-feIm Tr[P ri P rj Pr k ]. (45) 

This expression can assume nonzero values because the Cartesian components of the 
projected position operator PrP do not commute [31]. 

We have used (41o)-(41c) to evaluate the OMP contributions 5 LC , 5 IC , and a cs for 
finite samples, finding excellent agreement with the /c-space calculations using (346)- 
(34<f). As an example, the finite-sample results for the zz component are plotted as the 
symbols in figure 2. 



3. Linear-response expression for the OMP tensor 

In sections 2.2 and 2.3.1 expressions were given for evaluating M{£) under periodic 
boundary conditions. Used in conjunction with finite-field ab-initio methods for periodic 
insulators [36, 29], they allow to calculate the OMP tensor by finite differences. 
Alternatively, the electric field may be treated perturbatively [33]. With this approach 
in mind, we shall now take the £-field derivative in (1) analytically and obtain an 
expression for the OMP tensor which is amenable to density-functional perturbation- 
theory implementation [37]. It should be kept in mind that in the context of self- 
consistent-field (SCF) calculations the "zero- field" part of the Hamiltonian (4), 

^° = -!V + Vs CF (r-), (46) 

does depend on £ implicitly, through the charge density. As we will see, this dependence 
gives rise to additional local-field screening terms in the expression for the OMP. 

We shall only consider the case where the OMP is calculated for a reference state 
at zero field, which we indicate by a superscript "0." Upon inserting (34a) into (1) we 
obtain the three gauge- invariant OMP terms in (40). The term a cs is clearly of the 
isotropic form (2), with 



47ri 



9 CS = - — / d s ke ijk tr 



A°f) 4° - 2l A A A 



(47a) 



This is the same expression as obtained previously by heuristic methods [6, 7]§. The 
other two terms were not considered in the previous works. They are 

N 



LC 



ie abc f [dk\ J2 lm(2(d b u nk \(d c H° k )\d D u nk ) 

-(d b u nk \(d D H k )\d c u° nk ) (476) 

§ An inconsistency in the published literature regarding the numerical prefactor in (47a) has been 
resolved: see [38] 
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Figure 3. Contributions to the isotropic OMP from the Chern-Simons term a cs 
and from the Kubo-like terms a LC and 5 IC , expressed in terms of the dimensionless 
coupling strength 9 in (2). Model parameters are the same as for hgure 1. 



and 



522 = 7*-*/ [dfc] £ Im(2(« fe |9 D ^ fe )« fe |(9 c //°)K fe ) 

mn 

(47c) 



where <9£> denotes the field-derivative d/dSd and 

fewJifc) = <5fc I^Mnfc)| £=0 (48) 

are the first-order field-polarized states projected onto the unoccupied manifold. The 
terms containing doH^ describe the screening by local fields. They vanish for tight- 
binding models such as the one in this work, but should be included in self-consistent 
calculations, in the way described in [37]. We shall sometimes refer to 5 LC and 5 IC 
as 'Kubo' contributions because, unlike the Chern-Simons term, they involve first- 
order changes in the occupied orbitals and Hamiltonian, in a manner reminiscent of 
conventional linear-response theory. || 

Equations (47o)-(47c) are the main result of this section. The derivation of (476) 
and (47c) is somewhat laborious and is sketched in Appendix B. We emphasize that 
the Kubo-like terms, besides endowing the tensor a with off-diagonal elements, also 
generally contribute to its trace, which therefore is not purely geometric. Writing the 
isotropic part of the OMP response in the form (2), we then have 

e = e cs + e Kubo . (49) 

|| The terminology 'Kubo terms' for 5 LC and 5 IC is only meant to be suggestive. A Kubo-type linear- 
response calculation of the OMP should produce all three terms, including a cs . 
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Figure 4. Comparison between a zz calculated treating the two lowest bands as 
occupied (crosses) and the sum otzz + ol£z (thick solid line), where ai™^ (thin solid 
lines) correspond to treating only the lowest band (n = 1, upper line) or the second- 
lowest band (n = 2, lower line) as occupied. Model parameters are the same as for 
figure 1 except that the second lowest on-site energy in table Al is raised from —6.0 
to —5.0 in order to keep the two lowest bands well-separated. 



The two contributions are plotted for our model in figure 3. Moreover, the open circles in 
figure 1 show the zz and zy components of the OMP tensor computed from (47o)-(47c), 
confirming that the analytic field derivative of the magnetization was taken correctly. 

In the case of an insulator with a single valence band, the partition (40) of the 
OMP tensor acquires some interesting features. The terms a IC and a cs become purely 
itinerant, i.e., they vanish in the limit of a crystal composed of non-overlapping molecular 
units, with one electron per molecule. Also, the first term in the expression (47c) for a 
- the only term for tight-binding models - becomes traceless, as can be readily verified 
in a Hamiltonian gauge (where H k \u1 k ) = E^ k \u° nk )) with the help of the perturbation 
theory formula [33] 

\dnu° nk ) = ieJ2 'y^' lW- (50) 

m>N n m 

In order to verify these features numerically, we calculated the various contributions 
treating only the lowest band of our tight-binding model as occupied. The molecular 
limit was taken by setting to zero the hoppings between neighbouring eight-site cubic 
"molecules." 

It could have been anticipated from the outset that the Chern-Simons term (47a) 
could not be the entire expression for the OMP, based on the following argument [26]. 
Consider an insulator with N > 1 valence bands, all of which are isolated from one 
another. By looking at a da as dPd/dB a one can argue that, since each band carries a 
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certain amount of polarization , the total OMP should satisfy the relation 



where «W i s the OMP one would obtain by filling band n while keeping all other bands 
empty. We shall refer to this property as the "band-sum-consistency" of the OMP. It 
only holds exactly for models without charge self-consistency (see the analytic proof 
in Appendix C), but that suffices for the purpose of the argument. We note that the 
Chern-Simons contribution (47a) alone is not band-sum-consistent, because the second 
term therein vanishes for a single occupied band. Hence an additional contribution, 
also band-sum-inconsistent, must necessarily exist. Indeed, both 5 LC and 5 IC are band- 
sum-inconsistent, in such a way that the total OMP satisfies (51). This is illustrated in 
figure 4 for our tight-binding model. 

4. Summary and outlook 

In summary, we have developed a theoretical framework for calculating the frozen-ion 
orbital-magnetization response (OMP) to a static electric field. This development fills 
an important gap in the microscopic theory of the magnetoelectric effect, paving the 
way to first-principles calculations of the full response. While the OMP is often assumed 
to be small compared to the lattice-mediated and spin-magnetization parts of the ME 
response, there is no a priori reason why it should always be so. In fact, in strong Z 2 
topological insulators it is the only contribution that survives, and the predicted value is 
large compared to that of prototypical magnetoelectrics. Although the measurement of 
the 9 = 7r ME effect in topological insulators is challenging, as time-reversal symmetry 
must be broken to gap the surfaces [6, 7, 23], there may be other related materials 
where those symmetries are broken already in the bulk. The present formalism should 
be helpful in the ongoing computational search for such materials with a large and 
robust OMP. 

A key result of this work is a A;-space expression for the orbital magnetization of 
a periodic insulator under a finite electric field £ (equations (6a)-(6cf), or equivalently, 
(34a)-(34cf)). In addition to the terms (346)-(34c) already present at zero field [13], 
in three dimensions the field-dependent magnetization comprises an additional purely 
isotropic 'Chern-Simons' term, given by (34<f). This new term depends explicitly on £ 
and only implicitly on H%, while the converse is true for the other terms. Moreover, it is 
a multivalued quantity, with a quantum of arbitrariness M = £e 2 /hc along £. Thus, 
the analogy with the Berry-phase theory of electric polarization [10, 20], where a similar 
quantum arises, becomes even more profound at finite electric field. 

The Chern-Simons term M cs is responsible for the geometric part of the OMP 
response discussed in [6, 7] in connection with topological insulators. We have clarified 
that in materials with broken time-reversal and inversion symmetries in the bulk 
the CSOMP does not generally constitute the full response, as the remaining orbital 



N 




(51) 



n 
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magnetization terms, M LC and M lc , can also depend linearly on £. Their contribution 
to the OMP, given by (476)-(47c), is twofold: (i) to modify the isotropic coupling 
strength 6; and (ii) to introduce an anisotropic component of the response. 

Another noteworthy result is equation (45) for the Chern-Simons OMP of finite 
systems. One appealing feature of this expression is that it allows one to calculate the 
CSOMP without the need to choose a particular gauge. Instead, its /c-space counterpart, 
equation (47a), requires for its numerical evaluation a smoothly varying gauge for the 
Bloch states across the Brillouin zone. Equation (45) is also the more general of the 
two, as it can be applied to noncrystalline or otherwise disordered systems. 

We conclude by enumerating a few questions that are raised by the present work. 
Do the individual gauge-invariant OMP terms identified here in a one-electron picture 
remain meaningful for interacting systems, and can they be separated experimentally? 
(This appears to be the case for M LC and M IC at £ = [35].) How does the expression 
(47o)-(47c) for the linear OMP response change when the reference state is under a finite 
electric field El Finally, we note that equation (41c) for the CSOMP of finite systems 
has a striking resemblance to a formula given by Kitaev [39] for the 2-D Chern invariant 
characterizing the integer quantum Hall effect. Can this connection be made more 
precise, in view of the fact that the quantum of indeterminacy in 6> cs is associated with 
the possibility of changing the Chern invariant of the surface layers? These questions 
are left for future studies. 
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Appendix A. Tight-binding model and technical details 

Tight-binding model 

We have chosen for our tests a model of an ordinary (that is, non-topological) insulator. 
The prerequisites were the following. It should break both time-reversal and inversion 
symmetries, as the OMP tensor otherwise vanishes identically. It should be three- 
dimensional, as the geometric part of the response vanishes otherwise. Its symmetry 
should be sufficiently low to render all nine components of the OMP tensor nonzero. 
Finally, it should have multiple valence bands, for generality. 

We opted for a spinless model on a cubic lattice. It can be obtained starting 
with a one-site simple cubic model, doubling the cell in each direction, and assigning 
random on-site energies Ei and complex first-neighbour hoppings tj^i = te 1 ^^ of fixed 
magnitude t — 1. The Hamiltonian reads 




(A.l) 
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Table Al. The parameters of the tight-binding model. Columns I— III give the site 
coordinates i = (x, y, z), in units of the lattice constant a — 1 of the 2x2x2 primitive 
cubic cell. Column IV contains the on-site energies Ei, and the last three columns 
contain the phases of the complex nearest-neighbour hopping amplitudes along bonds 
in the negative x, y, and z directions. 



X 


y 


z 


Ei 


^(i + x/2) -> i 


+ y/2) -S- i 


4>{i + z/2) i 


0.0 


0.0 


0.0 


-6.5 


<pe [0,2tt] 


0.57T 


1.77T 


0.5 


0.0 


0.0 


0.9 


1.37T 


0.27T 


0.57T 


0.5 


0.5 


0.0 


1.4 


0.87T 


1.47T 


0.67T 


0.0 


0.5 


0.0 


1.2 


0.37T 


1.97T 


I.Ott 


0.0 


0.0 


0.5 


-6.0|| 


1.47T 


0.87T 


0.37T 


0.5 


0.0 


0.5 


1.5 


0.67T 


1.77T 


0.77T 


0.5 


0.5 


0.5 


0.8 


0.87T 


0.67T 


1.27T 


0.0 


0.5 


0.5 


1.2 


1.97T 


0.37T 


1.47T 



| In figure 4 the value —5.0 was used instead. 



where i = (x, y, z) labels the sites and (ij) denotes pairs of nearest-neighbour sites. The 
values of Ei on two of the eight sites were adjusted to ensure a finite gap everywhere in 
the Brillouin zone between the two lowest bands (chosen as the valence bands) and the 
remaining six. We also made sure that nonzero phases 4>j^i were not restricted to two- 
dimensional square-lattice planes, otherwise those are mirror symmetry planes, whose 
existence is sufficient to make the diagonal elements of the OMP tensor vanish. In our 
calculations all the model parameters were kept fixed except for one phase, which was 
scanned over the range [0, 2tt], and the results are plotted as a function of this phase ip. 
For reference, the on-site energies and the phases of the hopping amplitudes are listed 
in table Al. The energy bands are shown in figure Al for tp — 0. 

In order to couple the system to the electric field and to be able to define its orbital 
magnetization, the position operator r must be specified along with H° . We have chosen 
the simplest representation where r is diagonal in the tight-binding basis. 

Technical details 

The calculations employing periodic boundary conditions were carried out on an 
80 x 80 x 80 /c-point mesh, and the /c-space implementation of finite electric fields 
was done using the method discussed in section V of [30] . The open boundary condition 
calculations used cubic samples containing Lx L x L eight-site unit cells, that is, 2L + 1 
sites along each edge. For large L, we expect the magnetization to scale as 

M(L) = M + | + A + ^ (A.2) 

where a, b, and c account for face, edge, and corner corrections, respectively [13]. 
Calculations of M(L) under small fields were done using L = 4,5, 6, 7, and then fitted 
to (A.2) in order to extract the value M of the magnetization in the L — > oo limit. The 
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Figure Al. Band structure of the cubic-lattice tight-binding model given by (A.l), 
for the choice of parameters in table Al and ip = 0. 

differences between OMP values calculated in various ways as shown in figures 1 and 2 
were of the order of 1CT 7 e 2 / he or less. 

Before evaluating the /c-space expressions for M(S) [(Ga)-(Gd) and (346)-(34d)] 
and a [(47a)-(47c)] on a grid, they need to be properly discretized. The presence of the 
gauge- dependent Berry connection in (Get) demands the use of a "smooth gauge" for its 
evaluation, where the valence Bloch states given by (18) are smoothly varying functions 
of k. This is achieved by projecting a set of trial orbitals onto the set of occupied Bloch 
eigenstates according to the prescription in equations (62-64) of [31]. (For the tight- 
binding model discussed below, when treating the two lowest bands as occupied, the two 
trial orbitals are chosen as delta functions located at the two sites with lowest on-site 
energy.) If needed, this one-shot projection procedure can be improved upon by finding 
an optimally smooth gauge using methods based on minimizing the real-space spread 
of the WFs [31], but we found our results to change negligibly when performing this 
extra step. In a smooth gauge the needed /c-derivatives of the Bloch states and of the 
Berry connection matrix are then evaluated by straightforward numerical differentiation. 
Note that (Gb) and (6c) should be evaluated in the same smooth gauge as (Get), as these 
three equations are not separately gauge-invariant. A smooth gauge must also be used 
for (34c?) and (47a), because, as discussed in section 2.3.1, the Chern-Simons 3-form is 
locally gauge-dependent. 

The same strategy can be used to discretize (346) and (34c). However, since the 
k- derivatives appearing in those equations are covariant, the discretized form of the 
covariant derivative (19) given in [30, 13] may be used instead, circumventing the need 
to work in a smooth gauge. We have implemented both approaches, finding excellent 
agreement between them. 

Finally we come to equations (476) and (47c). In addition to the /c-derivative of 
the valence Bloch states, we need their (covariant) field-derivative (48), as well as the 
/c-derivative of H^. The latter quantity is easily calculated within the tight-binding 
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method, and for the former we used the linear- response expression (50). Note that this 
requires choosing the unperturbed states to be in the Hamiltonian gauge. This choice 
precludes calculating the /c-derivative on the right-hand-side of (50) by straightforward 
finite differences, which can only be done in a smooth gauge. But because (u Q mk \ddU^ k ) 
equals (Umkl^dUnk) f° r m > N, the discretized covariant derivative approach may be 
used instead. Alternatively, one can evaluate the ordinary ^-derivative by summation 
over states as 

\OdU nk ) - \ u mk) E o _ E o • l A ^J 

We note that this formula may not be used to calculate the geometric term (47a), 
because it induces locally a parallel transport gauge (A° n = 0), which cannot be enforced 
globally since the Brillouin zone is a closed space. 

Appendix B. Derivation of equations (47b) and (47c) 

For notational simplicity we drop the crystal momentum index k. So, for example, 
shall be denoted by \u n ). In order to calculate the OMP terms 

^ = d D Mi LC) \ £=Q (B.l) 



and 

532 = d D M^ 



(B.2) 



starting from (31) and (32), we shall first examine the field- and /c-derivatives of certain 
basic quantities. 

We begin by noting that the field-derivative dpP = —OdQ of the projection 
operator (20) can be written as 

N 

d D P = ^2 (\d D u n ){u n \ + \u n ){d D u n \} = d D P, (B.3) 

n 

in terms of the covariant field-derivative (48) (a similar expression holds for the k- 
derivative). This follows from a relation analogous to (21): 

N 

\d D u n ) = \d D u n ) - Aiu,d\ui), (B.4) 

i 

where 

Ain tD = i{ui\d D u n ) = A* nlD (B.5) 

is the Berry connection matrix along the parametric direction E&. With the help of (B.3) 
the field-derivative of (22) becomes 

d D F nmM = (d 2 Db u n \Q\d c u m ) + (d b u n \Q\d 2 Dc u m ) + i(F bD A c ) nm - i(A b F Dc ) nm , (B.6) 
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where F bD is obtained from F bd by replacing 84 with do- We shall also need the field- 
and /c-derivatives of the matrix H® m defined by (8): 

dn(H° nm )\ e = = 1 [ A d, H°] nm + (d D H° op ) nm \ £=Q (B.7) 

dc(H° m )\ £=Q = i [A° c ,H°] nm + (d c H° p ) nm \ £=0 , (B.8) 

where we introduced the notation (do jC H^ ) nm = (u n \dD iC H°\u m ) , where 'op' indicates 
that the derivative is taken on the operator itself, not its matrix representation. These 
two relations follow directly from (B.4) and (21). We will also make use of identities 
such as 

Re tr [XF bc ] = Re tr [X^F cb ] . (B.9) 
In particular, if X and Y are Hermitian, 

Re tr [XYF bc ] = Re tr [YXF cb ] . (B. 10) 

We are now ready to evaluate (B.2): 

5£ = -Itabc J [dfc]Im ti[F bc d D H° + H°d D F bc }\ £=0 . (B.ll) 

Inserting (B.6) and (B.7) on the right-hand side generates a number of terms. Some can 
be combined upon interchanging dummy indices b -B- c and invoking (B.10), leading to 

5^ = -l^abc J [dk] (2Retr [A D H°F bc + H°F bD A c ] + Im tr [F bc d D H° op ] 

N 

+21mJ2<n(d 2 D b u n \Q\d c u m )) . (B.12) 

mn 

Integrating the last term by parts in k b and using (B.3) and (B.8) again produces a 
number of terms, most of which cancel out. The end result reads 

Sffi = ie abc J [dk] Im tr [2F bD d c H° op - F bc d D H° op ] \ £=Q . (B.13) 

Similarly, (B.l) can be evaluated by repeatedly using (B.3) and integrating by parts the 
terms with mixed field- and A;- derivatives, yielding 

= Itabc J [dk] Im tr [2(d c H°) bD - (d D H°) bc \ \ £=Q , (B.14) 

where (d c H°) bD and (dDH°) bc are defined in analogy with (26), e.g., 

{d c H°)nrnk,bD = (d b u nk | (d c H%) \ d D u mk ) . (B.15) 

Equations (B.14) and (B.13) are respectively equivalent to (476) and (47c) in the main 
text. The gauge invariance of these equations follows from the fact that they are written 
as traces over gauge-covariant objects. (We also note that the covariant derivative 
transforms according to (18) regardless of the parameter with respect to which the 
differentiation is carried out.) 
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Appendix C. Band-sum consistency of the OMP 



Here we show analytically that the OMP tensor a satisfies the band-additivity relation 
(51) in models without charge self-consistency. In order to isolate the contribution 
coming from valence band n (assumed to be well-separated in energy from all 
other bands), we choose the Hamiltonian matrix to be diagonal at zero field, i.e., 
H^ nk (£ = 0) = E® k 5 mn . If in addition we use a parallel-transport gauge for the linear 
electric field perturbation [33] (this is achieved by setting to zero the matrix Ad defined 
in (B.5)) we find, using (B.7), 9oif^ nfe |f =0 = 0. With the help of these two relations, 
the field- derivative <9i)M a |£ = o of (6 a) is easily taken. From the first two terms therein 
we obtain (dropping the index k) 

N 

2-fe abc / [dk] ^(d b u n \d c (H° + E n )\d D u n ) 



(C.l) 



£=0 



In the parallel-transport gauge \dou n ) is given by (50), and combining the resulting 
expression with the field-derivative of the third term in (6 a) yields 



N 

a da = 2e 7 e a6c J [dk] ]T Re(d b u° n \d c (H° + E° n ) 



\1>N n 



!«?><«? 



f N 

e-fe abc / [dk] Re (( u m\ddU° n ) (d b u° n \d c u J) . 



(C.2) 



To find ctfa we replace Yi>n with Yi^ n and reduce sums Y mn anc ^ to single 



N 
mn 



terms. Inserting these expressions into (51) and splitting Yi^ n into Yi>n and X^n' 
some terms cancel and others can be combined, leading to 



^abc 



N N 

EE* 

n m^n 



X 



K\d d u° n ) 

(d b u:\d c (H° + £°)|<> 



0. 



(C.3) 



E° - E° 

^n ^rn 

The LHS is proportional to the difference between a da and Yin a da-> an d vanishes as a 
result of an exact cancellation between the terms (n, m) and (m, n) in the double sum. 
The integrand of the (n, m) term is 

' (d b u° n \d c (H° + E°)\u°J 1 



«\d d <) 



E° - E° 



+ ^(dbU° n \d c u. 



and after some manipulations the integrand of the (m, n) term becomes 



(u° m \d d u 



[<\d c (H° + E° m )\d b u° m ) 1 



The final step is to use the identity 
(d b u° n \d c (H° + E° n )\u° m ) _ (d b u° n \E° m -H°\d c u m 



E° - E° 

^n ^rn 



E° - E° 

^n 



d c (El + El 



(C.4) 



(C.5) 



(C.6) 
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(This identity follows from the relation 

(H° - E°J\d c u° m ) = -(d c H° - d c E°J\u° m ), (C.7) 

which in turn can be obtained by expanding H°\u^) = E^u^) to first order in the 
change in wavevector k.) The quantity (C.4) + (C.5) then becomes 

(d b u n \E m -H \d c ul) , (d c u° n \E° n - H°\d b u° m ) 



/..O 10 0\ / \yb^n\^m rc« m / , \^»nrn ^"W , /o Ola „,0 

\ u m\ u d u n/ \ fiO _ fiO ' fiO _ fiO ' \ u b a n\ u c a r. 



m 



.(C.8) 



Interchanging b •<-)■ c in the second term and combining with the first yields minus the 
third term, which concludes the proof. 
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